Load¶

Libraries and functions¶

Warning message in is.na(x[[i]]):
“is.na() applied to non-(list or vector) of type 'environment'”Warning message in rsqlite_fetch(res@ptr, n = n):
“Don't need to call dbFetch() for statements, only for queries”

==========================================================================
*
*  Package WGCNA 1.63 loaded.
*
*    Important note: It appears that your system supports multi-threading,
*    but it is not enabled within WGCNA in R. 
*    To allow multi-threading within WGCNA with all available cores, use 
*
*          allowWGCNAThreads()
*
*    within R. Use disableWGCNAThreads() to disable threading if necessary.
*    Alternatively, set the following environment variable on your system:
*
*          ALLOW_WGCNA_THREADS=<number_of_processors>
*
*    for example 
*
*          ALLOW_WGCNA_THREADS=4
*
*    To set the environment variable in linux bash shell, type 
*
*           export ALLOW_WGCNA_THREADS=4
*
*     before running R. Other operating systems or shells will
*     have a similar command to achieve the same aim.
*
==========================================================================


Allowing multi-threading with up to 4 threads.
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."

Data¶

Sample clustering¶

With all genes¶

Without X/Y genes¶

PCA on the normalized counts¶

With all factors¶

SPF vs GF¶

M, 8w¶

M, 52w¶

M, 104w¶

F, 8w¶

F, 52w¶

F, 104w¶

Gene co-expression analysis¶

Extract modules of co-expressed genes¶

pickSoftThreshold: will use block size 4347.
 pickSoftThreshold: calculating connectivity for given powers...
   ..working on genes 1 through 4347 of 10291
   ..working on genes 4348 through 8694 of 10291
   ..working on genes 8695 through 10291 of 10291
   Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k. max.k.
1      1    0.227  1.310          0.952 1030.000  1.03e+03 1570.0
2      2    0.232 -0.842          0.911  328.000  3.12e+02  758.0
3      3    0.590 -1.550          0.954  132.000  1.15e+02  446.0
4      4    0.729 -1.880          0.974   61.300  4.74e+01  292.0
5      5    0.771 -2.080          0.972   31.700  2.12e+01  205.0
6      6    0.812 -2.130          0.981   17.800  1.01e+01  150.0
7      7    0.802 -2.220          0.974   10.600  5.05e+00  114.0
8      8    0.831 -2.190          0.983    6.630  2.61e+00   88.7
9      9    0.852 -2.140          0.988    4.330  1.40e+00   70.5
10    10    0.867 -2.110          0.991    2.920  7.74e-01   56.8
11    12    0.872 -2.060          0.983    1.440  2.55e-01   38.4
12    14    0.877 -2.010          0.977    0.779  9.07e-02   26.9
13    16    0.893 -1.940          0.982    0.449  3.40e-02   19.4
14    18    0.906 -1.860          0.985    0.273  1.33e-02   14.2
15    20    0.929 -1.780          0.993    0.173  5.31e-03   10.7

Parameters for network construction

Soft-thresholding power parameter: 8 (R2 > 0.80 and mean connectivity highest as possible (Ref: https://www.degruyter.com/view/j/sagmb.2005.4.issue-1/sagmb.2005.4.1.1128/sagmb.2005.4.1.1128.pdf))
Minimum module size: 50
Height for module merging: 0.35 (below)
Biweight Midcorrelation

 Calculating module eigengenes block-wise from all genes
   Flagging genes and samples with too many missing values...
    ..step 1
 ..Working on block 1 .
    TOM calculation: adjacency..
    ..will use 4 parallel threads.
     Fraction of slow calculations: 0.000000
    ..connectivity..
    ..matrix multiplication (system BLAS)..
    ..normalization..
    ..done.
   ..saving TOM for block 1 into file norm_genes_TOM-block.1.RData
 ....clustering..
 ....detecting modules..
 ....calculating module eigengenes..
 ....checking kME in modules..
     ..removing 1 genes from module 1 because their KME is too low.
     ..removing 3 genes from module 9 because their KME is too low.
     ..removing 1 genes from module 12 because their KME is too low.
     ..removing 1 genes from module 13 because their KME is too low.
 ..merging modules that are too close..
     mergeCloseModules: Merging modules whose distance is less than 0.35
       Calculating new MEs...

 ME0  ME1  ME2  ME3  ME4  ME5  ME6  ME7  ME8  ME9 ME10 ME11 
4663 1239 1082  630  595  552  476  367  315  178  139   55

11 modules with between 55 and 1239 genes
4663 genes not assigned to a module

Relationship between modules and samples¶

Genes in modules¶

Percentages of XY genes per modules

[1] "0: 147 - 3.15247694617199%"
[1] "1: 36 - 2.90556900726392%"
[1] "2: 53 - 4.89833641404806%"
[1] "3: 18 - 2.85714285714286%"
[1] "4: 13 - 2.18487394957983%"
[1] "5: 19 - 3.44202898550725%"
[1] "6: 9 - 1.89075630252101%"
[1] "7: 7 - 1.90735694822888%"
[1] "8: 9 - 2.85714285714286%"
[1] "9: 2 - 1.12359550561798%"
[1] "10: 1 - 0.719424460431655%"
[1] "11: 1 - 1.81818181818182%"

Heatmaps with Z-scores¶

Type / gender / age¶

Gender / type / age¶

Sinaplots of the Z-scores per groups¶

The mean of the Z-score over the samples in the group is plot for each gene

Type / Age¶

Type / Age / Gender¶

Age / Type / Gender¶

Median of counts over ages¶

Enrichment analysis in modules¶

[1] "0"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
9622  GO:0034641            1.930318e-34                        1       1689
11864 GO:0044237            1.819262e-33                        1       2515
9636  GO:0034660            1.511191e-32                        1        212
11881 GO:0044260            2.808359e-32                        1       2109
2983  GO:0006139            6.603118e-32                        1       1518
9532  GO:0034470            7.640556e-32                        1        169
      numInCat                                             term ontology
9622      5285     cellular nitrogen compound metabolic process       BP
11864     8386                       cellular metabolic process       BP
9636       405                          ncRNA metabolic process       BP
11881     6844         cellular macromolecule metabolic process       BP
2983      4672 nucleobase-containing compound metabolic process       BP
9532       301                                 ncRNA processing       BP
[1] "1"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
17661 GO:0090304            3.630391e-20                        1        450
11881 GO:0044260            1.095788e-18                        1        644
5710  GO:0016070            1.118447e-17                        1        406
2983  GO:0006139            9.921890e-17                        1        468
12852 GO:0046483            4.487238e-16                        1        471
11364 GO:0043170            1.928903e-15                        1        667
      numInCat                                             term ontology
17661     4172                   nucleic acid metabolic process       BP
11881     6844         cellular macromolecule metabolic process       BP
5710      3729                            RNA metabolic process       BP
2983      4672 nucleobase-containing compound metabolic process       BP
12852     4765                    heterocycle metabolic process       BP
11364     7455                  macromolecule metabolic process       BP
[1] "2"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
12259 GO:0045184            9.868354e-13                        1        181
11879 GO:0044257            1.106652e-12                        1         84
14123 GO:0051234            2.182943e-12                        1        345
5312  GO:0015031            2.460062e-12                        1        167
13314 GO:0048193            2.934615e-12                        1         40
3431  GO:0006810            4.433883e-12                        1        333
      numInCat                                  term ontology
12259     1765 establishment of protein localization       BP
11879      604    cellular protein catabolic process       BP
14123     4062         establishment of localization       BP
5312      1606                     protein transport       BP
13314      185               Golgi vesicle transport       BP
3431      3915                             transport       BP
[1] "3"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
17661 GO:0090304            8.971982e-19                        1        245
5710  GO:0016070            1.893249e-17                        1        223
14136 GO:0051252            2.138832e-15                        1        190
2983  GO:0006139            5.935538e-15                        1        251
12852 GO:0046483            7.973750e-15                        1        254
14083 GO:0051171            1.403879e-14                        1        210
      numInCat                                              term ontology
17661     4172                    nucleic acid metabolic process       BP
5710      3729                             RNA metabolic process       BP
14136     3046               regulation of RNA metabolic process       BP
2983      4672  nucleobase-containing compound metabolic process       BP
12852     4765                     heterocycle metabolic process       BP
14083     3601 regulation of nitrogen compound metabolic process       BP
[1] "4"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

Network visualization¶

TOM calculation: adjacency..
..will use 4 parallel threads.
 Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

Power	SFT.R.sq	slope	truncated.R.sq	mean.k.	median.k.	max.k.
1	0.2273633	1.3069408	0.9522773	1032.2103062	1.031743e+03	1574.40467
2	0.2318240	-0.8419065	0.9109768	328.3750082	3.124206e+02	758.22758
3	0.5898524	-1.5468341	0.9541521	131.7106717	1.147058e+02	446.45478
4	0.7293365	-1.8760203	0.9740987	61.3157190	4.736691e+01	292.33277
5	0.7706350	-2.0775848	0.9724795	31.7270126	2.117472e+01	204.76919
6	0.8116126	-2.1271367	0.9807292	17.7731786	1.008260e+01	150.23963
7	0.8022879	-2.2224796	0.9735875	10.5925115	5.050971e+00	114.00564
8	0.8307511	-2.1879600	0.9827800	6.6345171	2.608895e+00	88.74112
9	0.8519550	-2.1423593	0.9884310	4.3282040	1.397797e+00	70.46121
10	0.8665930	-2.1052573	0.9914710	2.9212252	7.738615e-01	56.84323
12	0.8724663	-2.0611029	0.9830544	1.4447620	2.552006e-01	38.38129
14	0.8766131	-2.0101797	0.9767692	0.7788890	9.066889e-02	26.89092
16	0.8933757	-1.9356365	0.9817629	0.4488075	3.402853e-02	19.36386
18	0.9058564	-1.8638594	0.9845667	0.2726832	1.328237e-02	14.24265
20	0.9285263	-1.7822036	0.9931630	0.1729982	5.314736e-03	10.65525

	0610007P14Rik	0610009B22Rik	0610009O20Rik	0610010K14Rik	0610030E20Rik	0610037L13Rik	Card19	1110008F13Rik	⋯	Zswim7	Zswim8	Zufsp	Zwint	Zxda	Zyg11b	Zyx	Zzef1	Zzz3	Hikeshi
0610007P14Rik	1	1	0	0	0	0	0	1	⋯	0	1	0	0	0	0	0	0	0	0
0610009B22Rik	1	1	1	0	0	0	0	0	⋯	0	1	0	0	0	0	0	0	0	1
0610009O20Rik	0	1	1	0	0	0	1	0	⋯	0	1	0	0	1	0	0	0	0	0
0610010K14Rik	0	0	0	1	0	0	1	0	⋯	1	0	0	0	0	1	1	0	0	1
0610030E20Rik	0	0	0	0	1	1	0	1	⋯	0	1	1	1	0	1	0	1	1	0
0610037L13Rik	0	0	0	0	1	1	0	0	⋯	0	1	1	0	0	0	0	1	1	0